############################################################
# Unpacking urban (dis)continuities of postwar violence    #
#                      Replication code                    #
############################################################


## Load packages needed for replication (install if necessary)

library(lme4)  
library(lmerTest)
library(performance)
library(stargazer)
library(tidyverse)


## Open the dataset (update text to your local directory as appropriate)
load("C:/..../replicationdata.Rdata")



## REPLICATION OF MAIN MODELS (Table 1 in article) #######


# model 1
mlm1 <- lmer(
  GED_log ~ 
    warintensity_city_5year_re + 
    pop_log + 
    urbrate + 
    econ_importance_dummy + 
    PKO_presence_city + 
    war_country_SB_re + 
    GED_fatalities_country_re + 
    secondary_warring_5yr_binary + 
    PKO_presence_country + 
    (1 | country), 
  data=replicationdata)

summary(mlm1)
r2(mlm1) # inspect R2 for mixed model - conditional R2 reported in article
ranef(mlm1) # inspect country-level intercepts


# model 2
mlm2 <- lmer(
  GED_log ~ 
    national_capital + 
    regional_capital +
    warintensity_city_5year_re + 
    war_country_SB_re +
    GED_fatalities_country_re +
    pop_log + 
    urbrate +
    econ_importance_dummy +
    secondary_warring_5yr_binary +
    PKO_presence_country + 
    PKO_presence_city + 
    (1 | country), 
  data=replicationdata)

summary(mlm2)
r2(mlm2)

# model 3
mlm3 <- lmer(
  GED_log ~ 
    national_capital+
    regional_capital+
    incomp_govt +
    national_capital*incomp_govt +
    warintensity_city_5year_re + 
    war_country_SB_re +
    GED_fatalities_country_re +
    pop_log + 
    urbrate +
    econ_importance_dummy +
    secondary_warring_5yr_binary +
    PKO_presence_country + 
    PKO_presence_city + 
    (1 | country),
  data=replicationdata)

summary(mlm3)
r2(mlm3)


# model 4

mlm4 <- lmer(
  GED_log ~ 
    national_capital+
    regional_capital+
    incomp_terr +
    regional_capital*incomp_terr +
    warintensity_city_5year_re + 
    war_country_SB_re +
    GED_fatalities_country_re +
    pop_log + 
    urbrate +
    econ_importance_dummy +
    secondary_warring_5yr_binary +
    PKO_presence_country + 
    PKO_presence_city + 
    (1 | country),
  data=replicationdata)

summary(mlm4)
r2(mlm4)

# model 5

mlm5 <- lmer(
  GED_log ~ 
    victory_reb+
    victory_gov+
    negotiated+
    warintensity_city_5year_re + 
    war_country_SB_re +
    GED_fatalities_country_re +
    pop_log + 
    urbrate +
    econ_importance_dummy +
    secondary_warring_5yr_binary +
    PKO_presence_country + 
    PKO_presence_city + 
    (1 | country),
  data=replicationdata)

summary(mlm5)
r2(mlm5)


class(mlm1) <- "lmerMod"
class(mlm2) <- "lmerMod"
class(mlm3) <- "lmerMod"
class(mlm4) <- "lmerMod"
class(mlm5) <- "lmerMod"


# Write model output to a regression table using Stargazer
stargazer(
  mlm1, 
  mlm2, 
  mlm3, 
  mlm4,
  mlm5,
  title="Table 1: Postwar urban violence (GED)", 
  dep.var.labels= c("postwar violence (GED)"),
  notes.append = FALSE,
  star.cutoffs = c(.1, .05, .01),
  omit.stat=c("rsq","f"),
  type="html",
  out="results_mlm.html"
)





## REPLICATION OF APPENDICES   ###################

## Appendix B: Descriptive statistics 

df_descriptive <- replicationdata %>%
  select(
    GED_total_fatalities_city, 
    warintensity_city_5year,
    warintensity_country_5year,
    GED_fatalities_country,
    urbrate, 
    population,
    econ_importance_dummy, 
    national_capital, 
    regional_capital, 
    victory_gov, 
    victory_reb,
    negotiated,
    PKO_presence_country, 
    PKO_presence_city,
    secondary_warring_1yr_binary,
    secondary_warring_5yr_binary
  )

stargazer(as.data.frame(df_descriptive), type="html", Title="Descriptive statistics", digits=2, out="summarystats.html")
stargazer(as.data.frame(df_descriptive), type="text", Title="Descriptive statistics", digits=2, out="summarystats.txt")




##  Appendix C: termination, different reference category 

mlm5a <- lmer(
  GED_log ~ 
    victory_reb+
    victory_gov+
    lowactivity+
    warintensity_city_5year_re + 
    war_country_SB_re +
    GED_fatalities_country_re +
    pop_log + 
    urbrate +
    econ_importance_dummy +
    secondary_warring_5yr_binary +
    PKO_presence_country + 
    PKO_presence_city + 
    (1 | country),
  data=replicationdata)

summary(mlm5a)
r2(mlm5a)

mlm5b <- lmer(
  GED_log ~ 
    victory_reb+
    negotiated+
    lowactivity+
    warintensity_city_5year_re + 
    war_country_SB_re +
    GED_fatalities_country_re +
    pop_log + 
    urbrate +
    econ_importance_dummy +
    secondary_warring_5yr_binary +
    PKO_presence_country + 
    PKO_presence_city + 
    (1 | country),
  data=replicationdata)

summary(mlm5b)
r2(mlm5b)

mlm5c <- lmer(
  GED_log ~ 
    victory_gov+
    negotiated+
    lowactivity+
    warintensity_city_5year_re + 
    war_country_SB_re +
    GED_fatalities_country_re +
    pop_log + 
    urbrate +
    econ_importance_dummy +
    secondary_warring_5yr_binary +
    PKO_presence_country + 
    PKO_presence_city + 
    (1 | country),
  data=replicationdata)

summary(mlm5c)
r2(mlm5c)

class(mlm5a) <- "lmerMod"
class(mlm5b) <- "lmerMod"
class(mlm5c) <- "lmerMod"


stargazer(
  mlm5a,
  mlm5b,
  mlm5c,
  title="Robustness: Postwar urban violence (GED)", 
  dep.var.labels= c("postwar violence (GED)"),
  notes.append = FALSE,
  star.cutoffs = c(.1, .05, .01),
  omit.stat=c("rsq","f"),
  type="html",
  out="AppendixC.html"
)



